close all; clear all; path(pathdef); clc;

% add subfunctions
addpath(genpath('subfunctions/'))
addpath(genpath('data_figure1/'))

% load cell divistion times measured as a function of global temperature
% it takes ~10 seconds to finish
% load_cell_division_excel_data;
load('worm_division_temp.mat')

% color library
color_list{1} = [67 84 147]/255;        % dark blue
color_list{2} = [110 153 201]/255;      % light blue
color_list{3} = [196 40 27]/255;        % dark red
color_list{4} = [218 134 121]/255;      % light red
color_list{5} = [40 127 70]/255;        % dark green
color_list{6} = [161 196 139]/255;      % light green
color_list{7} = [218 165 32]/255;       % gold
color_list{8} = [238 221 130]/255;      % light gold
color_list{9} = [128 0 128]/255;        % purple
color_list{10} = [204 153 255]/255;     % light purple
color_list{11} = [161 165 162]/255;     % grey


%% Perform clustering of data %%%

% threshold for group difference screening
thrsh=0.15; 

cell_measT_cluster = cell(1,14);
cell_measT_cluster_error = cell(1,14);
cell_division_time_cluster = cell(1,14);
cell_division_time_cluster_error = cell(1,14);

for jj = 1:14
    [datax,idx] = sort(cell_measT{jj});
    datax=datax(~isnan(datax))';
    datay = (cell_division_time{jj}(idx));
    datay = datay(1:size(datax));%%Exclude Nan cases corresponding to non-measured cases
    nGrps=sum(abs(diff(datax))>thrsh)+1;           % number of groups with that difference or more
    C=clusterdata(datax(:),'maxclust',nGrps);      % do the cluster with those groups C is index to grp
    cnts=accumarray(C,1);                         % the counts of each group
    
    cell_measT_cluster{jj}=accumarray(C,datax(:),[],@mean);          % compute means of each group
    cell_division_time_cluster{jj}=accumarray(C,datay(:),[],@mean);
    cell_measT_cluster_error{jj}=accumarray(C,datax(:),[],@std)./sqrt(cnts);  % compute error of each group
    cell_division_time_cluster_error{jj}=accumarray(C,datay(:),[],@std)./sqrt(cnts);
end

%% Arrhenius law - two cell %%%

bPlotOption.MarkerSize  = 6;
bPlotOption.CapSize     = 0;
bPlotOption.LineWidth   = 1;
bPlotOption.FontSize    = 12;
bPlotOption.MarkerList  = {'o', 's'};
marker_idx              = [1 2] ;
color_idx               = [5 1];

% which cell?
cell_idx_list = [1 2];

close(figure(200));
fh = figure(200); cla()
set(fh, 'position', [100 100 [230 200]])

for zz = 1:length(cell_idx_list)
    
    % cell selection
    cell_idx = cell_idx_list(zz);
        
    % data processing
    xdata = cell_measT_cluster{cell_idx}+273;
    ydata = cell_division_time_cluster{cell_idx};
    ydataerr = cell_division_time_cluster_error{cell_idx};
    ydataerr(ydataerr == 0) = min(ydataerr(ydataerr>0));
    
    temp = [xdata ydata ydataerr];
    temp = sortrows(temp);
    
    xdata = temp(:,1);
    ydata = temp(:,2);
    ydataerr = temp(:,3);
    
    % fitting 
    fitrange = 2:length(xdata);
    FitType = fittype('1e-15*b*exp(1e3*a./x) + 1e12*d*exp(-1e3*c./x)');
    FitRes = fit(xdata(fitrange), ydata(fitrange), FitType, 'startpoint', [10 1 1 1], 'lower', [0 0 0 0], 'weights', ydataerr(fitrange), 'robust', 'BiSquare');
        
    % plotting
    xfit = linspace(10, 30, 100) + 273;
    plot(xfit, FitRes(xfit), '-', 'linewidth', bPlotOption.LineWidth ,'color', color_list{color_idx(zz)}, 'handlevisibility', 'off')
    hold on
    plotcustom(xdata, ydata, ydataerr, bPlotOption, marker_idx(zz), color_idx(zz))
    
    bReading = [283 288 293 298 303];
    set(gca, 'yscale', 'log')
    set(gca, 'XTick', bReading)
    set(gca, 'XTickLabel', bReading-273)
    set(gca, 'YTick', [20 40 60])
    axisinfo('Temperature (degC)', 'Cell cycle time (min)', [12 28]+273, [10 60], bPlotOption.FontSize)
    
end

legend boxoff
legend('AB','P1')


%% Arrhenius law - four cell %%%

bPlotOption.MarkerSize  = 5;
bPlotOption.CapSize     = 0;
bPlotOption.LineWidth   = 1;
bPlotOption.FontSize    = 11;
bPlotOption.MarkerList  = {'o', 's', 'd', '<'};
marker_idx              = [1 2 3] ;
color_idx               = [5 1 3];

% which cell?
cell_idx_list = [6 5 4];

close(figure(201));
fh = figure(201); cla()
set(fh, 'position', [100 100 [150 200]])

for zz = 1:length(cell_idx_list)
    
    % cell selection
    cell_idx = cell_idx_list(zz);
        
    % data processing
    xdata = cell_measT_cluster{cell_idx}+273;
    ydata = cell_division_time_cluster{cell_idx};
    ydataerr = cell_division_time_cluster_error{cell_idx};
    ydataerr(ydataerr == 0) = min(ydataerr(ydataerr>0));
    
    temp = [xdata ydata ydataerr];
    temp = sortrows(temp);
    
    xdata = temp(:,1);
    ydata = temp(:,2);
    ydataerr = temp(:,3);
    
    % fitting 
    fitrange = 2:length(xdata);
    FitType = fittype('1e-15*b*exp(1e3*a./x) + 1e12*d*exp(-1e3*c./x)');
    FitRes = fit(xdata(fitrange), ydata(fitrange), FitType, 'startpoint', [10 1 1 1], 'lower', [0 0 0 0], 'weights', ydataerr(fitrange), 'robust', 'BiSquare');
        
    % plotting
    xfit = linspace(10, 30, 100) + 273;
    plot(xfit, FitRes(xfit), '-', 'linewidth', bPlotOption.LineWidth ,'color', color_list{color_idx(zz)}, 'handlevisibility', 'off')
    hold on
    plotcustom(xdata, ydata, ydataerr, bPlotOption, marker_idx(zz), color_idx(zz))
    
    bReading = [283 288 293 298 303];
    set(gca, 'yscale', 'log')
    set(gca, 'XTick', bReading)
    set(gca, 'XTickLabel', bReading-273)
    set(gca, 'YTick', [10 20 40 60 80])
    axisinfo('Temperature (degC)', 'Cell cycle time (min)', [12 28]+273, [10 90], bPlotOption.FontSize)
    
end

legend boxoff
legend('P2', 'EMS', 'ABa/ABp')


%% Arrhenius law - eight cell %%%

bPlotOption.MarkerSize  = 5;
bPlotOption.CapSize     = 0;
bPlotOption.LineWidth   = 1;
bPlotOption.FontSize    = 11;
bPlotOption.MarkerList  = {'o', 's', 'd', '<'};
marker_idx              = [1 2 3 4] ;
color_idx               = [9 7 5];

% cell types

% which cell?
cell_idx_list = [12 11 10];
% cell_idx_list = [14 13];

close(figure(202));
fh = figure(202); cla()
set(fh, 'position', [100 100 [150 200]])

for zz = 1:length(cell_idx_list)
    
    % cell selection
    cell_idx = cell_idx_list(zz);
        
    % data processing
    xdata = cell_measT_cluster{cell_idx}+273;
    ydata = cell_division_time_cluster{cell_idx};
    ydataerr = cell_division_time_cluster_error{cell_idx};
    ydataerr(ydataerr == 0) = min(ydataerr(ydataerr>0));
    
    temp = [xdata ydata ydataerr];
    temp = sortrows(temp);
    
    xdata = temp(:,1);
    ydata = temp(:,2);
    ydataerr = temp(:,3);
    
    % fitting 
    fitrange = 2:length(xdata);
    FitType = fittype('1e-15*b*exp(1e3*a./x) + 1e12*d*exp(-1e3*c./x)');
    FitRes = fit(xdata(fitrange), ydata(fitrange), FitType, 'startpoint', [10 1 1 1], 'lower', [0 0 0 0], 'weights', ydataerr(fitrange), 'robust', 'BiSquare');
        
    % plotting
    xfit = linspace(10, 30, 100) + 273;
    plot(xfit, FitRes(xfit), '-', 'linewidth', bPlotOption.LineWidth ,'color', color_list{color_idx(zz)}, 'handlevisibility', 'off')
    hold on
    plotcustom(xdata, ydata, ydataerr, bPlotOption, marker_idx(zz), color_idx(zz))
    
    bReading = [283 288 293 298 303];
    set(gca, 'yscale', 'log')
    set(gca, 'XTick', bReading)
    set(gca, 'XTickLabel', bReading-273)
    set(gca, 'YTick', [10 20 40 60 80])
    axisinfo('Temperature (degC)', 'Cell cycle time (min)', [12 28]+273, [10 90], bPlotOption.FontSize)
    
end

legend boxoff
legend('E', 'MS', 'AB_{al,ar,pl,pr}')